{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Large-Scale Stochastic Variational GP Regression in 1D (w/ KISS-GP)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This example shows how to perform GP regression, but using **variational inference** rather than exact inference. There are a few cases where variational inference may be prefereable:\n",
    "\n",
    "1) If you have lots of data, and want to perform **stochastic optimization**\n",
    "\n",
    "2) If you have a model where you want to use other variational distributions\n",
    "\n",
    "KISS-GP with SVI was introduced in:\n",
    "https://papers.nips.cc/paper/6426-stochastic-variational-deep-kernel-learning.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gpleiss/anaconda3/envs/gpytorch/lib/python3.7/site-packages/matplotlib/__init__.py:999: UserWarning: Duplicate key in file \"/home/gpleiss/.dotfiles/matplotlib/matplotlibrc\", line #57\n",
      "  (fname, cnt))\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a training set\n",
    "# We're going to learn a sine function\n",
    "train_x = torch.linspace(0, 1, 1000)\n",
    "train_y = torch.sin(train_x * (4 * math.pi)) + torch.randn(train_x.size()) * 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing SGD - the dataloader\n",
    "\n",
    "Because we want to do stochastic optimization, we have to put the dataset in a pytorch **DataLoader**.\n",
    "This creates easy minibatches of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "train_dataset = TensorDataset(train_x, train_y)\n",
    "train_loader = DataLoader(train_dataset, batch_size=64, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "This is pretty similar to a normal regression model, except now we're using a `gpytorch.models.GridInducingVariationalGP` instead of a `gpytorch.models.ExactGP`.\n",
    "\n",
    "Any of the variational models would work. We're using the `GridInducingVariationalGP` because we have many data points, but only 1 dimensional data.\n",
    "\n",
    "Similar to exact regression, we use a `GaussianLikelihood`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.GridInducingVariationalGP):\n",
    "    def __init__(self):\n",
    "        super(GPRegressionModel, self).__init__(grid_size=20, grid_bounds=[(-0.05, 1.05)])\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.RBFKernel(\n",
    "                lengthscale_prior=gpytorch.priors.SmoothedBoxPrior(\n",
    "                    math.exp(-3), math.exp(6), sigma=0.1, transform=torch.exp\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "    \n",
    "model = GPRegressionModel()\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The training loop\n",
    "\n",
    "This training loop will use **stochastic optimization** rather than batch optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/20 - Loss: 1.485\n",
      "Epoch 2/20 - Loss: 1.314\n",
      "Epoch 3/20 - Loss: 1.203\n",
      "Epoch 4/20 - Loss: 1.070\n",
      "Epoch 5/20 - Loss: 1.028\n",
      "Epoch 6/20 - Loss: 1.030\n",
      "Epoch 7/20 - Loss: 0.752\n",
      "Epoch 8/20 - Loss: 0.728\n",
      "Epoch 9/20 - Loss: 0.479\n",
      "Epoch 10/20 - Loss: 0.355\n",
      "Epoch 11/20 - Loss: 0.216\n",
      "Epoch 12/20 - Loss: 0.124\n",
      "Epoch 13/20 - Loss: 0.135\n",
      "Epoch 14/20 - Loss: 0.068\n",
      "Epoch 15/20 - Loss: -0.013\n",
      "Epoch 16/20 - Loss: -0.001\n",
      "Epoch 17/20 - Loss: 0.053\n",
      "Epoch 18/20 - Loss: 0.015\n",
      "Epoch 19/20 - Loss: -0.059\n",
      "Epoch 20/20 - Loss: -0.004\n",
      "CPU times: user 19.6 s, sys: 852 ms, total: 20.5 s\n",
      "Wall time: 5.11 s\n"
     ]
    }
   ],
   "source": [
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},\n",
    "    {'params': likelihood.parameters()},\n",
    "], lr=0.01)\n",
    "\n",
    "# Our loss object\n",
    "# We're using the VariationalMarginalLogLikelihood object\n",
    "mll = gpytorch.mlls.VariationalMarginalLogLikelihood(likelihood, model, num_data=train_y.numel())\n",
    "\n",
    "# The training loop\n",
    "def train(n_epochs=20):\n",
    "    # We use a Learning rate scheduler from PyTorch to lower the learning rate during optimization\n",
    "    # We're going to drop the learning rate by 1/10 after 3/4 of training\n",
    "    # This helps the model converge to a minimum\n",
    "    scheduler = torch.optim.lr_scheduler.MultiStepLR(optimizer, milestones=[0.75 * n_epochs], gamma=0.1)\n",
    "    \n",
    "    for i in range(n_epochs):\n",
    "        scheduler.step()\n",
    "        \n",
    "        # Within each iteration, we will go over each minibatch of data\n",
    "        for x_batch, y_batch in train_loader:\n",
    "            optimizer.zero_grad()\n",
    "            output = model(x_batch)\n",
    "            loss = -mll(output, y_batch)\n",
    "            \n",
    "            # The actual optimization step\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            \n",
    "        print('Epoch %d/%d - Loss: %.3f' % (i + 1, n_epochs, loss.item()))\n",
    "            \n",
    "%time train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f38401cff98>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADDCAYAAABtec/IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJztnXl4VOX1+D+zZs9MJjt7EvY1gVFBQYWACFitCG5Vab8o7c9ataW4ILjUfaFWrdIiqLhQLRFtNYgssYoalIGA7FvCTtbJZLLPdn9/TCYGmC2ZSTIh7+d5eEjm3rn35J57z33f855FJkkSAoGgeyLvbAEEAkHnIQyAQNCNEQZAIOjGCAMgEHRjhAEQCLoxykAPoNfrJzf9OMVgMDwY6PEEAkHHEdAIoOnhn20wGDYCo/V6/ejgiCUQCDoCWbDiAPR6/RGDwZARlIMJBIIOISg+AL1e/wDw22AcSyAQdBzBHAGsBu4yGAwmd9sfeughEXIoEHQSzz33nMzd5wE5AV1zfoPBsB0oBOYBL3ja/4knnvB5zNLSUpKSkgIRq90JdRlDXT4IfRlDXT7wX8bHHnvM47ZApwCTAV3Tz1qcRkAgEHQRAjUAy4B0vV4/D8BgMOQELpJAIOgoApoCNM33lwVJFkE3wWazUV1dTXV1NaGajepwODCbzZ0thlfOlVEmkxEWFkZKSgpKpX+PdsCBQAJBaykuLkaj0RAfH49M5tY31elYrVZUKlVni+GVc2WUJAmTyURxcTG9evXy6xgiFFjQ4TQ2NhIbGxuyD39XRSaTodVqaWxs9Ps7wgAIOhxJkkLi4S8oKKCgoKDdz2MymVizZk27nwecRqA10yphAAQhzZkzZ5g8eTLFxcVtPkZBQQHLly9n06ZNLF++nMJC52KVRqMhJ6f9/dZardbteQoKChgyZAhr1qxhzZo1LFmypFk2d3jb1laED0AQ0jz77LN8//33PPPMM7z66qut/r7JZOLFF19k1apVzZ/deuutrFq1Cp1O5+WbwSUuLu68z7KyskhLS2PmzJnNn02fPp21a9eet29hYSErVqzg6aefDqpcwgAIQhKtVktDQ0Pz78uWLWPZsmWEh4djMrkNNnVLTk4OkyZNOuuzuLg4Nm3axJgxYygoKGDTpk3s2LGDuXPnsm3bNgC2bdvGrFmzyMvLQ6fTkZaWRlFRETk5OaSlpTFo0CDWrVvHqlWr+P3vf8/8+fMBzto/LS2NFStWkJmZyfbt2/3+u11v+ry8PAAmTZrEjh07KCoqoqCgAI1GQ15eHna7nSlTppCenu739TgXMQUQhCT79u3jpptuIiIiAoCIiAhuvvlm9u/f3+pjVVVVedyWlZVFdnY2mZmZrFixgh07dpCXl8fEiRNZtGgRY8aMaX74J02aRFxcHE8//TR33HFH8zFmzpxJenr6efs/8sgjXH/99WRnZ5OWltYqmdPT09HpdOh0Oj755BMmTZpEWloaWVlZ520LBGEABCFJamoqsbGxNDY2Eh4e3rxykJKS0qrjTJo0qfmt7qKoqIjs7OyzPnNNB66//nrmzp3LkiVLsFgsaDQasrKymkcRWq32rGMvWbKEMWPGNH927v6txWQykZ6ezpIlS9BoNGRmZjZ/Ds6pgGvbqFGjztrWFsQUQBCylJaWctdddzF37lxWrFjRJkdgeno6CxYsYPny5aSlpbFjxw7+/ve/N283mUxnTQFcQ/aJEycyZcoUVqxY0fz2dQ3BTSYTWq2WWbNm8cgjjzQbhaeeeuqs/efPn88nn3xCZmZm83ezsrKaz11QUEBRUVHzCkFRUVGzbK7zVVVVUVhYSGVlJSaTiaKiouZtRqORwsJCioqKzjpuawhaNqAvHnroIUkkA3UMoS7f4cOH6du3b0gH2nTFQCAXhw8fpn///s2/P/bYYx6zAcUUQCDoxggDIBB0Y4QBEAi6McIACATdGGEABIJujDAAAkE3RhgAgaAbIwyA4IKmoKCAcePGnZX2W1hYeN5n3ZVgtAab1/RjhmgNJmgt4eFhQTlOQ4P7IhhZWVnNkYCvv/464MwNcMXVd3eC0Rpso8FgcBUHnezrOwJBR6PRaDxuKywsZPny5axZs4aCgoLm39966y0KCwvZtGkT06dPZ9OmTTzyyCMdKHXHEOgUIB1naXBwlgRve16ioFvS0NAYlH++mDlzJsuXLz8vHv/cDL5zM+2ys7PRarVkZ2cHlHQTqgRaFbhlReDRwEeBiSMQtA/Z2dnceuutZ2XuudBoNKSnp5OWlsaSJUvIzMykd+/eHD9+HJPJ5LaYx4VCULIBmzoEbW/qEOSR0tJSn8fqClY21GUMdfkcDgd2u71DzrVjxw7efPNNevfuTWZmJr169WLr1q0UFBSwdetWHn/8cZYtW8akSZPo168fffr04fDhw5SXl3P48GE+++wzCgsLOXjwIIWFhWzdurU5Rbez8XQNHQ6HX88aBCkbUK/XP2AwGDy2BAORDdiRhLp8IhswOIRENqBer5/neviFE1Ag6FoEYxXgeb1ef0Sv11cGSSaBQNBBBOoE3AhcuB4SgeACR0QCCgTdGGEABIJujDAAAkEb6Mh2X+2JqAos6HRe++pIQN//w8QMr9sLCgrYtm1bc23+vLw8vzvsuAKDduzY0dz8A35u99Wyq09XRBgAwQWNu9Zg/r65TSYTRqOR7Oxst23ELoQIQWEABBc0OTk554X/zp8/n8LCQvLy8khLS6OqqgqNRsOSJUuYP38+eXl5PP7442zbto2ioiI2bdrEokWLyM/Px2Qyndfuy3UsV0swo9F41rGefvrp5t5+LlkyMzPP+k5nZSYKH4Cg2+Fq43XnnXeSnZ1NTk6O26QfV5JQdnY2o0ePBnDb7uvchCJ3x1qyZAlz585l5syZTJo06bzvdBbdbgRgqrNibrBSZ7FT02hHkiSGpMQQHd7tLkWnYrVLWCUbNrtEvcWOTCZDLgOFXIZcJkPmNnC19cyaNYu77777rM82bdoE0Nzhx2QyBZz00zKhCNxPD1zTCFcnoXO/0xl0m7v+RGU9245VcqKy/rxtPx6tZFBKNKN7a9FFqTtBuu5DaXUjtRY7VQ025HLnANQuSeDKSbGDDBlKhQy1Qh6wIdBqtWe1BquqqiIzM5OnnnqKnJwcdDodM2fOpKioiKKiouZWWzt27MBsNje3Atu+fTsFBQVu232d2xLs3GO5vvfiiy82v/XP/U7LnoMdyQXfGuxYRR0/Hq2k2Nzgc18ZMgalRDNpUCIKuec7L9STbUJVvp9OVvHtkQrGxVtI6dWn2QB4Qi6TEalWEKlWIAvWkMBPuksy0AU7ArDYHHx9qJz9xdV+f0dCYn9xNQ1WB9OGJaFUCBdJMLDYHGw6UMbh0ppWfc8hSdQ02qi32tGEq1AphT6CzQV5RU+bGli19WSrHv6WHK2oZe2eEmx2R5Al635IksQXe0pa/fC3xO6QqKy30mDtmBoC3YkLzgBsKTKypuA01Q3WgI5zrKKOz3eXYBVGICDyC40cN9ad9ZmE0zC0BkmSMNfbqLPYgijdhYckSa2aLl1QBiC/0MjWo5VIBMevccJYx9rdJTgcHeMnudA4VFrDtuPnVyeqtcmor61uvRFAorrBRk2DrdXf7Q5IkoTJZCIszP9KyxeMD2D7cROGY8EvSXDcWEd+kZHLMuKDfuwLmYoaC5v2l7ndtsekwGY3ojGZaKtrL0wpR92OPgGHw+HTSdnZnCujTCYjLCyMlJQUv49xQRiAvWfMfHekot2Ov/24iaSYMAYkRbfbOS4kGq12cncXe5w+WSUZP5yxERUV1eZzyJAxdVhSu+kkVFdSWhIMGUPbxPnBkbJa8g6Ut/t5Nu0vo6LG0u7nuRDIO1hOVX1gPhhfSEhs3FdGidl3SXCBZ7q0ASirbuTLvSUdMh+02h3k7i6mUXiivXKwpCYgj39rsDkcrN1dTE2jcAy2laAYgKay4B1Ko9XOF3tKsHegg66q3soGD/NaAdQ02vjfwfYfjZ17zrW7xJJtWwlGVeDJwOogyNIqNh0oa/Uw02ws440FczAb2/4QF5XXsqe4zveO3ZC8A2U02lo/QgpULyXVDXx3xNim73Z3AjYATYVBC4Mgi98UnDBxpKzW537n3lgbVy3l6J5tbFy11O12fzGcrKa8Rsw9W7LntJljFf4ZRrOxjLcf/W1Q9fLTqSoKy33fE4Kz6XKrAKdN/lt714311O2TkByuIWJP8nOPkJ/7AJAKjCN3+WpueeBuL0c6G7tD4ss9pdyo74lKhAtjrrfy7WH/V2HWv7ecY3sVPHnbeyANB+4EZpGf+w35ubcAe5DJHGxctZSZ9zzq93E37S8j6aIwosO63G3daXTolQq0NVitxc5ne43UWbwPM5+6ZTw2awuPvSQB04FHgEvP23/7Vza2f/UvFMqXuX/p0+S8/Aiz/vg0MXEJbo/f0NDAibJKPt/WyGX9Yn3+TR1NR7cG+2K/kUqz7xUSp16uA/4JxHF+vNaspv+NSNIa8nMfJz93GAqlil4Dh3vVCUAt8PEPh7l6UFzAyUOh3l4NgiNjhxoAf9cs3e1nsTnYuP00MlU4UR6StMzGMt5/9s/c8/K/+DrnLXZ9n4fNMgNkj4Dk6udWBRxBFWbGZj2B5IgAfgncjt12OysWbsdUVs33n670+vaJioriZC1Uy6LISGz7enZ70VFr2LtPmzHbVUR5UgpOvbz79EKGXLKPXd+6GkjvRybfSURUEXXV3yFX9MJhvxS5YhIOe0+co4JbSe77X1L7fcrOb/7tUycAZjucaFCj7xt4ua5QjwOAwGXsEmMlu0Mid3cxFbXe592uIf8PX/wbiX7YLF8CE0ACmbyEPoM2cM2dMWzP+5D9hs2YSk+jVKuxWVKBPwFzMZWNBr4hP3cS+bnDUKrUPPvfAo/n/OpAGb204YSpFMH8k7sENQ02vvNj6P/Zsk85tvd1IB2Z3Irk+BMK5TIcdithkT0YdfkExk6fwpa1q9m39WFMpTHI5H9Bcsym5NjNlBybCISTn7uS/NyPfOpkS1ElveMiSY71PyS2uxKwAdDr9bOc/+lnGQyGnCDIdB55B8o46aaQh4uHr806a8ifn6sEngY0RMY2kNznAyKiP+Y3j70EQL+hi1n55H0M1k9g7PTZbFm7GmPJatTqfHbn34UkTQS+ZvBFT3D1nGzeWDCH2x5+iVhd4nnnrrfa2VJUyRUDPQ9NL1TyDpZh8bL85tTLMOBrIAY4jOS4EZlsJ3c99y47v/qc6spyZt6zGICZ9zj1MuSiBMZOV7D27Qc5vOMG7LaLgXeQK0YzcsIGrpz1G686kSSJTfvLuEnf02tdB0GIFwSx2h0YjvmO8Tcby/j8zRfZ9b0Bm+Vl4BYABunLuOXPRqI0/i1NffzaE2xZ+zky2adI0hSUKjNDL3maXd+9yNjpNzYPP2tra88KY5XJZNw0pieJMaHxxumIMNb9xdVs2Ofdp3PmqInX7u+PtTEJ+ASl+reMuOxirrlzAYqwSL9CgXNefYIfvogGlgJqdCk/0D/z72z98oOzdOKOsWk6LurXtqnAhRQK3KUKguQdMmE/2khto93vNeVYXSIKVSI2yxdAFlBD+ogPmPuXCa0qKVVjMjJuxrXop9Ty2v252Kwz+OnbB4EN5Od+1Dz8XPSvb8/6niRJfH2onBuyenR45ZrOoLbRxjeHvA/97Xb4zz9GYW2MAn5Aofo1dms14ZHRxOoSqa31b8mutsrIuBlyemT8wMevjsRYfAk/rrMBX5ylE3dTgq3HKslIjBJl3rwQcgbgTLUFZZh/82mX0+/GP/6VvVvmA8PRJlXTb+hibJY9yGQTWnXuOYtfaf554cpyXp+fT1X5OOA/KNWXMeKykVxz5wL3clc1sL+4hiGpMa06Z1fk2yMVXo2z2VjGK/eexFwxFKXKyMjLV3LFzLfZsnY11ZWtixRsqZOEHgd567F+WBsvAzajVGcz4rLRHnVid0h8dbCcmZmp3cIwt4WQMwCtYeOqpRTt/ol/PBRDXfVwNAlW7n7xDHFJ8wI+dlxSAoP0T/HjOgdwGTbLSsLC3/T69vruSAXpCZEXtEPwTFUDB0u8x/r/64VdmCvuQSazMe+ZKtKG3wPQPNdvK/1HRTHissVsz/s9MAKb5V3U4W+79QO4OG2qZ/fpakb0DL3l2lCgS0axPHxtFgumDSM/dzXwPlVlmUApNaYxxCUFLzGkzlzGmOwPiIxtACZyaOf150WxtaTeamfL0eDXJAgVJElis5ehv1MvN3J459ym/f/AGwvSefja4DW9sDQcYfSkfxIZ0wBM5lDBbJ+Rg98fqaCmQSQMuaNrGoB31pN15Qzk8teB2UAVg/RPsXDla0E9z5zFr3Dzn+/h14vPIFdIVJyezYcv7eD4vh3NYavnsuuU+YJNGz5QUkNJtefqyg8s30BE9KdABPA2SvU7ZE2cwcPvrA+aDHMWv8ItC+7mN0+cQaF0YCz+Je8/e+isUOJzsdgd5BeJXAF3dEkDEKtLxFw5FYfjd0ADcA265FKvQ8FASBteD9IfAThUcCeSNIz83I9YMG3YeW83SZLYfLhjM+I6AovNwfc+QrB//HII9TX9gUIUqgXYrY3NTr9g029IPZJ0JwBFu+cgSVd41AnAgeIaSqtF/sa5dCkD4BrqHf6plsJdTuVn37KXcTN6UV3ZfhWBABa+ez265E1AFPBvlGqdx7fbicr6Cy4xZdtxE7UeCnKajWW8fM9zbPowAXAwdOz73Pu3ZYydflO76uWRd+eQ1DsHUAE5KNX9PepEQmrXqlFdlS7lBHQ6/fbz7l96IDkiGJ1tYurt4chkgTmX/EETn0hG5rMYv0wBhmGzPEt4ZJ7Ht9u3hyvoq4u8IAJRzPVWCk54jjtf/95yTh95GJAx4Xoj185zxvQH6vTzRawukbThSyk9EQlMx2Z5k7CIf3rUycnKeorKa0lLCL3Q7c6iS4wAfnb6fQQso762J7CbnV+nBa2HnD/UVxczcsK7KJR2YB4nDw/xuG9VvZUdJ6s6Trh25LtCo9vCKy69/LBuJDAU2M/mT3oF1enni9qqCvSTP21y1F5J4e6rvO7//RGjqPLcgq5hAFxOP+V9OKP8ahhy8d9ZuPLTDpVjzuJXuH3hHK664yQA5acWceKQ2aMH2nC0ktouXq7qRGW9xxJfD7+zngFZfwbmA3YUqrvImpgdVKefL+YsfoWb5t/Prx4sAaD85K/Z+0OjR50Y6yzsOdO2hjGhhLE2OI7mLmEAYnWJ1NeOw2F7qemTu9AmmtrkXFLK5cSGq0iOCadffBRpCVHIWzmMuHhaKUMurqa+RsE7j8dStLvArQe6q3ufHQ7vy34RUUmcOPhHQI5M/hIO23d+O/0iVQr6xkei7xvH1cOSuWJgAv2ToolUt21WOnB0LVfcUI7DIWPV8z0p2n3I46rAj0crsdi6bgmxBqudtbtLgnKsLuEDOFMUxsHt9wNK9FMOoFLTaudSTLiKMX20DEmJPq/nX3WDjR0nq9hz2uxXJyCZDA5uHwwYMBuHAgvIz33ObVjq/jM1DEmJoac2olXyhgK7T5u9ZmCufz+RhtoEIqJPcudTgzBsuMlnpF+kWslVQxIJs9WcF8c+sqcGcL7dfjplZvdpc6sKvn77n37A1zTWjwFeJz/3Nrc6qbPY2HbcxLh0nd/HDhUcDmerNVO9FQg89yTkDUBVhZK3HuuDw65i1OVVzL7fjlzuv3NJG6FC3zeOQcnRyD045GLClUzoH8/FfbVsO25y283mXBauXMWq59/gyE9PAo+iVH/KiMsyzgtLlZDI21/GLRf16lLNRustzixHTxw/EM7Xa+KRySXufNJKn0ED6TPIu156x0UwZUgSUWFKSr1UDtZFqblyYAIje8by7ZEKv0uNLVyZS84r/2Tfj38FfoVC+TkjJ1S7DRUuOGFiSEoM2sjQ7gB8Lt8cruBkZX3QQptD+o4sP13BC3faMZWp6Du0jpvmn8bfZi0yZIzuo+XWi3szJDXG48PfkjCVgksz4rl8gO/U3lhdIom9DgLvAxHYLM97HP6a6q382MUiBLcUGT3G+xtLyvnHgxKSQ8bl11fQZ7DnVG1wZkuOTdNx3ahUolpRrksXpebakan8clQqMX58L1aXiCahCngAALvtZZSqZLc6sTskvu1iy4K7TpnZdSq4juWQGQE4HLB0qYK1m3tiNkZQVa7i1JEeWBu1hEUW85tHq1Cp/RsORocpmTwkid5xbRt2j+qlQS6T8fXBcq99BmtMRsZkf8eu727C0nAtp45s9bjv9hNV9E+KJilEUoa9UVbdyJ7Tnh1l7z5VhrXxcsKjzjD1dt+jpUmDEhia2vZY/N66SGaN6clnPxX7LMZaYzIydrqCEwfKOXUkhUM7Znvct6i8lmMVdfSNj2yzbB3FmaoGvjkU/ACzkDEAcjn85S8Kqqpa9jWLBI7TWDeZx28+5LMSDEBGYhSTBiUSHmBCzoiesSjlMjbtL/NoBFyZaj0zKvjvshRqTIuxNh5BFXb+/pLknArcOKanX6ORzsLhkPjqgHvD5yzw4fR9gIOG2tks/OV3XvUyoqcmoIffRXSYkhuyerBub4nXKYFLJ6UnTbx8tw5T6dUcLDjGwCz3gVmbD1fQKy4ipOM1LDYHX+4txdEOtTtCagpw3312pt52iuv+3z76Zz6EUj0UyEAVdsKvmPJByTFMG5Yc8MPvYkhqDNmDfXu0L73WSEq/BozFata9G+FxCaqsppHtJ0I7NuDHY5Ue4/0fens90Zr3cUbevYEqbJtXvfTQRjChf/CaqqqVcq4ZnsKwHr4NSlIvC5Nvdb4xV7+cxN/nu0/gqqyz8FOQh9XB5quDZQG3u/dESBmAhQvtXDmrmPHXSiT0KMJu3Y9SrcBm8R1TnpYQxeTBiUHP+x6SGsPIXhqv+ygUcP3dxQBs/rQHRbuNXpagjEFbww02xVUNGI55HtIf3ZtOTdUIoAyF6kmveokOUzJtWHLQ36xyuYxJgxL9MgJXzionNb0BU1kEx/Ze71EnW4+afFaa7iz2nan2mX4dCCFlAFrinMvdxB9eXuUzprx3XARXD01qt6H1+Ix4UjXhXvd585FBwHtIDjXwN4+JKXaHs15dqPW3t9gcbNhX6lEua6OM3OXJAKQP/w/3/m2pR70o5DKmD08hUt1+dRGuGJDgc2l10cwszhSOB+zAfeTn7nerk0abnf8dDL2Wb5V1Fr5uh3l/S4LRGmyWXq+frNfrHwj0WGfOnGHZwrswG8uYs/gVZt6zmB7pg5l5z+KzKsO0JCU2nOnDU9p1iU0hlzFtWDJRXoJUHn5nPcPHfQmYgV+gUF3jcXhcbG6gIMSmAt8dqWhaWz4fs7GMF+b9j8pSNalpDfz2uUu96mV8Rny7V+R16SQ23PMynjOCNAW54g1AgUz2MplXutfJkbJadp0yt6PErcPukFi/t9RjXIq5oozZs2dTXFwc0HkCempcTUGb2oOZAm0S+uyzz3LUS679ucSEKblmRApqZfsPZKLClFw9LNlj1GCsLpHoOCvwDAB26/OowzUepy1bikJnKnC0oo7dpz3f/LlvfYip9A4Arv1tMXIvL/ZecREdVn0nQq3gFyNTUHsw/rG6RMIio3DYH8fZbGQiDbXjvSZwhUothx+PVnpNX96waik//vgjzzzzTEDnCfTJuQlwTRoLgcltOYhWqyU8PJxly5YhSZLXvG4XcpmMacOTiWjHYea59NCGM96LU6vGZOSSq0+hSagBhnJ8/ziP+7qmAp2dmFLTYGOjh+q+rmSf7ZsuxZkG/TH/fCjNo17UCjnZg4Lvh/GGLkrNVUOTkOH+nM5Cr1dxxQ3O/I2i3bdg9zDdtzkcrNtb4lc0aHtSYm70GIz2c2Lch0iSxLJlywgPD0er1bbpXIEuA2qBlsHuXl2+nlqDbd68maeeeoovv/yShoYGlOowhlwykavuuNdj/b1L+sQiazBT2tCxw7ZUNcSr7VS4kWvWn5zWuN/IM3z0wgAqy35HWfFPRMa4v+OO1NaSp7QwMjW46an+toyyOyTW7jdSXuN+6H/fG5/yyWtfUvjT7UAjCtUiho692qNeMvvF0lBdSYMfuTbBbL0VBWRoYOdpzzqxWS389G0DlSV9+fa/R9Ff5X7OX1tby+fbGhmm7RwjYHdI/HdPBTX17pPI7nvjU75c+QoHfvwfVksj4eHhXH311SxatMiv1nvnEhKtwZKSkkhKSsJisaBUqbFbLUTHakjp1dft/v2Topk4LLk9RfXK1GEONh5t9FggY8yVVrZ9WcvhnVHkrUqm4swvPTaxOGiSGJWuCXpPAX/qxX99sJw6SU2Uh7LZERFRlBxzFvSUyV/FYTtAdOwot3rpo4tkwvDUoMvoL1MTEqmRTlNs9lyy7Jq5Zbz3TG+++rAHuzbP4Y5Ff3Grk5O10CM2nIGd0Bfg+yMVWORhREW5vx+ioqKIjtVgtaShVtdhsZwkKSmJ4cOHt+l8gU4BTIAro0ILtDm2srS0lLvuuou7X3jHq9dfG6Eie1D7lP7yl3CVnClDEj0OO2Uy51xZJpcwbEyhaHe9R7+G3SGxbk9Jh2enHSip8bn+vW2Thtqq/qjCjNz94kCPeglTKjpdJ3K5jKuGJnn0BwCMGF9N3yF11JrVHNt3lVdf0+ZCc4eXECsxN/gVJ1JZIkMdvpm4uH3cfPNDlJS0PTMw0BHAR4C+6ed0YGObD/TRRwC8vG63x0oyKoWcacOTO8Tp54veukiy+mjY7mGu9up9Q5EcLwN3A6+QnzvZYxMLU72VvANlXN1Bo5qKGgtfHfC+7NVQK2ft2055bri3nn5DB9BvqHu9TOgfT3R45weVaiJUXDEwwWPHooXXZWGzjgbygT+RnzvQY/9Hm8NB7q5iZo/p2SHtxm12Bxv3+V4erq+VU1XxDpaGcHoPb+T11x8lIoBE04CeJIPBsB1Ar9dPBkyu39sDGTKmDk0iITp0YunHpuk8xvY//M56Rlz6Lc5BUTYK5Z1eo+YOlfp+IweDmkYbubuLPS8vNdVdXPt2FDWVSvoOqWP0RM9y9dFFhlQzlMEpMQxMjna7zbksGI9MngNEIJO5+oQXAAAVPUlEQVS/5FUnNY02cnd5vlbBQpIkNh0ow1jneQXCbCzj9fl38tajyRQfDSepdyMrVpQH9PBDEOIADAbDMoPBsNFgMCwL9FjeGN8/PuRquSnkMqYMSXK7NBirSyRKawec1YTttueRK3p6jWb89nBFuw476y12/rPzDFUe1vvBVXexmi1rU5HJJH75/4o9ll1TK+RM6uShvzsmDkwkxk18gGtZUHL8GahHctyEpWGMV52UVjeywY83cyBsKar0Ge234YN/cHTvHzi6N46YOCt3PnmCuLjADVPnj6X9YERPDZm9vYfjdha6KDWj+7hfgnFGM1rpO6QE0HGo4Favx3L5A9ojLLXR6nz4PcUenF138a9IkhJJWsHr8z3XPbw0I56YEBj6n4taKWfSIPcp3c5lwXFcMu0oAEd+ug2Hj8t9pKyGzYcr2sUI7D5t9tr81qWXLWtH4yyHZ6a68iJemNc2p9+5hLwB6KOL5PIgJpS0B/q+WrdvnDmLX+GGPyzm1gfNqMMdmCuuYOt6h9cuNlX1Vj7efgqzl7d0a7HYHHy2q5gyL6m0rrqLCuVMYDpQxfDLvvI4PO4VF8HwHqEz9D+XPrpIBqecL58rwvTaeaBJsNJQO4DNnyq86gRg58kqNuwrc1scta0crajjfwe9h/o+/M56+g19BlgAWFGobiFrYi8WBqnuYsgZgH66cMb00XLlwER+MTKVacOSQzp9FpzOSW9Zb7pkK1fPcTqmPnmjh9d6deB0CuZsPx2UqLSaRhuf7yrmTJXn5TFwDo/lymTstr83ffI4MVq72+GxqmnoH+oNNyf0jyfCQ2aoOlzimrlO7/m6d1Mo2n3YZwTqgZJqPttVHJQVm1OmetbtKfE5qrBZe3Di4L0AyBULcNi+CGqzlZAzAJf1i+XSjHhG9IylX3xkSHj8/SEjMcprYYncFWnAFqyNOuAFn9GOtRYbHxec5rTJ+4Prjf3F1az68SSnTN4r9gBIEhzaPgdIpUdGOWOnl3hcir24XxyaiNAvpRWuUnDFQM/VnT5ckgF8g82iARb7FYF6wljHmoLTba723Gi189WBMj4pOOPTuWi3warne2K3RRGXvJV7X7k86M1WQm8C14W5vH8CqypPuB0mLly5jtV/e5f9W0cD81Co/sPI8TKPra3BmaX2n51n0PfVMrxHrN9hz3UWO5/vKqbIj+5ErhbrI8e/i9l4OWERduYsqkSXssjt/toIFaN8pEeHEgOSojlYUuO2U9PClev5919zOLBtPPAHlKoPGDE+1atOwFnX4V9bT5LZW8PInhq/X1KHS2v45lCFxwCylpiNZbx6/xGqyoaiSbBy36uxRMUODnqzla7xeu0iaCNVHh2CsbpEtIkVwBMA2K1LUaoSfQ7lbA4HW4qMvJN/nK8OlFHpYanIYnNwqLSGDftK+XhXhV8PP7i8/mV8vtwZ3Xfd74rRpXj2P4zvHx/S1XPcccXABMKU5xvPWF0iccklwDJAhc36NqqwOL+G1/VWO/mFTr1sKTJS78FxW1bdiOFYJau3neKLPSV+PfwAH7/6DVVlvwYc3PrAKaJi26degRgBBJkxfbTsP1NNtZshYo3JyNhpRynaU0nJ8T4c9lKv7lxsDge7T5vZc7oaTYSSMKUCtUqOWiHHYrNzytTQXDLKanfgPrj3Z5zlvSyADNiE3RYJfMLHr93MRVe5L+/VNz4y5JZi/SE6TIm+r9Ztb8Aak5GLp27hYMGvMJWO5PCOWa06dqPNztajlRiOmghXyYlQKYhQKwhTyimtbqSmlVMFp16igV0438+Ps/SBJ/wqh9cWxAggyKgUci7NcO8QnLP4FW649xFuX1iBQumgsnQ627+y+fRAt0RCwlRvpaS6gRPGOo6U1XCisr7V9eKauy0pHgYmAiUMv2w1C1e69y7LZTLGe/i7ugKjemnc1g6Ys/gVZt//AHMWl6FQShiLr2PregdvP+q+hJgnJCTqrXaMdRZOmZzNYVv78INTL3FJHwM9gM0o1d6DlQJFGIB2YGBytNcKQsl9LUy9w3lzffxqT7880MEmVpdIQ/1FOOx/afpkLjFah8fh74ieseg8JA11BRRyGZdmeG4E0qt/A9P/z7kqsObvfTm2t6zDdQLOsmuVpVcCtShUv8VurWu3FusgDEC7MaF/gsdkIYAv3+0H5GNpiAde9ssDHUyqypUc3H4voOCiqw4wbkaMR+9yhErBJf3iOkSu9mRAUjQpsZ4N89q3+gK52KwxwHvk567uUJ3UmBSsec2ZUZk27CPu/dtT7d5iXfgA2onk2DAGp0Szr9h9cvzClevIeeUd9v2YCfwahfJrRk4o8+mBDhSzsYz3nn4Im3UddquG/pk1zLrXjlzh2bs8Nl1HWJAqLXc24/vHk7P9lNttC1eu55PXl7H7+9HAlcgVLzLq8rx21wlAVUUZf73bRp15EP1H1XLXM2ORy9u/xboYAbQjY9N1qLyUq9IkmID7AbDbXkNyDGy3oZ6LjauWcnTvLZw8pEGTYOVXD57yWuIrMTqMYSGU7BMoqZpwMhLdJwvF6hKJ1tqAWwErDvufqCr/RbvrBOCjl3ZSZ85Grqjnxj/63wErUIQBaEdc3mdPOHMFzAzWHwei2bf1D1ScqWiVU9Bffo71VwL3Ahaqysfz9B0jvX7v8gEJIR/x11ouy9B5rO3ozBVI4YrZ2wAo3D2XnZut7aITcOllOod2/AoAh/1envn1QJ/TjrggBWIJA9DOZPbSoIt07zxz5QrctrCWxF6NNNalseJRG0f3bAu6A+rhd9aTPvwJ4G0A5MoFZE1M9OpdHpAUTQ+t93LoXRGNl2AmV67AxJsUXDGrHCQl/3ohg6LdNe3iFLz/9TzCI7/BWU0vF6X6fZ9ef02EiusyW1d9yRPCALQzSoWciYM8Vw8CCIuQMBZfCtRTdnIiknRb0JyCrvz+ot1aivY+BKiQyV/EYXvNq3dZKZd36WU/X1zcL85jnoCLbz/tBXyC3RYNfEZ+7qagOQVd+f2rXx5GQ10f4CcUqt9gt3pvghMTruL6zB5BK1IiDEAH0EMbznAfpbIXrvwrfQa5SiosQ6G8MSjrv85IPxWrXhiE5FCT3Hc9972awbgZ3r3LY/pqQ6LKT3uhVsoZl+55WRBg4covGTnhfZBtB/qD7DuGjv2/oKzJb/hgKUf33s2xfVqUqgpGT/on9/7tH169/tFhSq7PTA1qCvaFq+EQ49J0HUfLa91GCILTAZWa/j3HD8iAe7HbVmGueBOANxbM8VhU1BM/R/rpgQ047OHASspP/ZaeGdu9epdjwlWMDtH6C8FkaGoMu06ZPaZJx+oSiYxRgjQDZOtAGsWhgpeorjQCJ3n/2T8HoJcngduAamzWyfy0eS+3LPi9R72EqxRcNyo16ElYYgTQQaiVcq+ZaQC1VUbGTs9n7PR9gIIjP/2OlU+ebJVPwGws4+1Hf8v/e+lDUvu9B3wHaJHJPyXzio9ZuPJLn8e4YkB8u3ZaChVkMhmXD/CuE6dT8ErufvEYsfF7sDbGsXRBPz5+9ZtW6+WNBXP49WOfoEn4DlgE2FAobyNrYk+vowoZMiYPTmyXQKygaDnQjkDdhbSEKAYkuV+Cgp+dgjf8QUIuvx9wcHz/r5CkV5vnnw/9IvMsj7TrxnL9vnHVUo7ttfLW4lGcOXoboEYm/weS42YioiN9vq2GpsZ2yXj/ttJDG+5TJzPvWUzasP48/JYcmWw1jfUK9v7wIJL0OPm5G/zSy4YP/kHR7rG8/fiVVJVfClQjV9yJw/6Zz0i/Ub017aaTYPQGnAysDoIs3YLLByT4dD4BPPLe9fQd8hJgBe4BjhOfuo5BY+ad9ebZuGopR/ds48nbZrBg2gLyc68EDNSa04GjwGTufy2JcTOu9xlRFhuuCmo7767CZRnxKP1YeFeqJR55N4nEnv8B1MCjwDGSeucwbOwt5+mlaPc2nrztdyyYtowta+cDrzUnXcFw7nv1Wp+Rfskx4Vzqw1cRCAH7AAwGw0a9Xl8YDGG6A5FqBVcOTOCLPd5rucfqEklN28axfZcjky1CkmZQcWYqFWemAv9Hfm4Z+bl7cHZjux2YBPxcIlYuf5Phl63nut8tJlaX6DOiTCZrqqvfRQqwBJOYcCWXpMW5zRY8F01CIhmj3qDs1EtNeplK6YkbKD0xA3iQ/Nxy8nP3AzOAZ0FKa/HtkyiUf2LkhDquufNdn3oJUyqYOiypXdOvO9QJ6E/romC2jGovApUxFkgKs1Nk9F7tx1Rein5qT/RTTHz177so3DkOa+NNgPsZl0K5E7vtU+TKz5Hs21BHzEQRFumxvVpLRqZGoWg0U1raMa3WQk3PPcMktEobp6qcDsGGBs+6ceolGf2Ucjauuo+ju6Zis04HhrrZuwx1+GYsDZ+gUObisJtQqPzTy8UZWhqrKyn10GotGNcwJFqDtXW/ziRQGa+Li+eDH09S56VAxNwnXm/+OeOJLD5+7Qm2rP0TCuUo7LZINPEDqKpoRK6w4rCvIzZezmD9BEZNvJ+dX31OdWU5UVG+546J0WFMzerZ4YU+Qk3PM7XxfGg41awTT9fuLL087dLLHSiUfbDbNGgShlFV3ohCuQ+7bQtRmhTGZE9g7PTlbFm72i+9DEmN4ZLBvq9PoNfQpwHQ6/Xz3Hxc2NQSXNBGwlUKJg5KIHeX//3dnR7paYydPpsta1ezO/9Dxs3Ibvq9nurKcmbes5ja2loy/EwiUSvkXDW0fYeZXYWoMCVThiTy353+6wRcepnSQi//btLLn5ofeNdQ35/knpgwJZf39746ESxkwah1rtfrNxgMhine9nnooYekJ554wuexSktLQ+7NcC7BlHHDvlL2e8gYbCu1tbV+vfkBpg9PISOx473+oaznbw9X8O3+U35fw2Bz3ahU+ug8F5h14e81fOyxx3juuefcWvhgrALMcv6nb10tJQEAl/ePd9tToCMY3UfbKQ9/qDMuXUdCVOfoZHiPWL8e/mARjFWAHCAnCLJ0S8JUCq4emsTHBadbXdYrEHrFRTAurf2Wl7oyCrmM7AFa/nfc0qayXm0lNlzFZR2cf9H91nxCkBRNuNdyVcEmJkzJ1KGh33ClM4lSK5gxIsVjPYdgI0NG9uDEDl+GFQYgRMjqrSW9AyLw1Ao504anEOlnj4HuTFJMGFcNTfKayRksLk6Lo1dcgK1+24AwACHE5MHuu9oGC6VczowRKSTHhk6L9VAnPSGK8e0cHTkkJYaLO6nmojAAIUSYSsG0dor8kstkTBue3Clvma5OZm8Nmb09V3YKhF5xEUzsxBbrwgCEGMmx4Vw7MhV1EOeerjDffl56Fwq8M6F/PGOD7DTVRamZPiy5U2MwhAEIQXrFRXB9Vg+/koZ8oZA7U0m9ZbwJ/OOifnFB64ocHabkFyNSOr3asjAAIUpSTBg3jO4RkE8gOSacm/W9GJxy4VT17WyG9Yhl+vBkv7IH3SFDxrAesdx6US9iQ6DDsqgIFMLERaqZldWD/x0q52h5HRL+xQnIZTLGpesY3VsrlvragfSEKGaP6cEPRZUUtUIv8VFhTByU4LVrVEcjDECIEx2u5JoRKZjqrPx0qop9Z6qxeOgrr4tU0zc+kkRlBIP6dv1OPqFMQnQYM0akUF7TyNZjJo6U1ro1BGFKBb3iwukXH8ng5JiQM8jCAHQRtJEqLh+QwNg0HaerGrDYHDTaHFjsDsKUcvrqIpuLRfqTdi0IDgnRYUwblow53Yq5wUaD1U6DzYHV7iA5JpyU2LCQe+hbIgxAF0OtlAtvfggSG6EKiTl9axFOQIGgGyMMgEDQjREGQCDoxggDIBB0Y4QBEAi6McIACATdGGEABIJujDAAAkE3JuBAoBZlwzMMBsODgR5PIBB0HAGNAJr6Am40GAzLgPSm3wUCQRch0ClAOs7mdACFTb8LBIIuQkBTgKY3v4vRwEfe9he9ATuGUJcPQl/GUJcPQqg3oF6vHw1sNxgM273tJ3oDdhyhLh+EvoyhLh+ETm/AycIBKBB0PXwagHOG+eeh1+vnGQyGF5p+niyahgoEXYdgrAI8r9frj+j1+sogySQQCDqIQJ2AGwFRe0og6KKISECBoBsjDIBA0I0RBkAg6MYIAyAQdGOEARAIujHCAAgE3RhhAASCbowwAAJBN0YYAIGgGyMMgEDQjREGQCDoxggDIBB0Y4QBEAi6McIACATdGGEABIJujDAAAkE3RhgAgaAbIwyAQNCNCUZrMFdjkCmiMrBA0LUIRlHQ2U21AUc39QcQCARdhGAUBXWVAU/31RhEIBCEFsHqDPQA8Ftf+z322GPBOJ1AIAgSMkmSgnIgvV6/GrjLYDCEflM1gUAABNgazDXnbxr6FwLzgBeCK6JAIGgvAm0NNhlwzfu1wNZgCCUQCDqGgKYAer1eC9zY9OsYg8Hg0w8gEAhCh6D5AASdg16vnwWYgNGuJq0e9nvA23ZB6KPX60d7Wmnz9z44l6CsArQVX0K39Y/qQPlc/pGMzgiCauGD2ajX69M93SBN8RpT6AT/jB/XcDSQDmAwGHI6WDyXDP7eh+m+umW3F006/CeQ4WabX/eBOzotFLil0IDp3CAiX9tDQL7JwMamGyK9RURkR3ITzhsTnE7YzpDBI37q8OGmBz+9MwLJ/LwPC5u2F3ZWsJvr/B42t/k+6MxcAF9Cd/bN7ev86S0+K2z6vaPRAsYWv8efu0PT22DjuZ93EF6vYdObdSuAwWB4oZMCyfy5z55v+j9Ug9183gee6EwD4EvoNv9RQcLr+Q0Gw7IWw8HRgKGjBGsluk48ty8dXgTE6/X60U3BZJ2BLz1vx/nmrzxnvwsCkQ0YIE1Dwu2d9GYw8fMDrgUqWm7s5Le/v1S4rl3TiCCkaFrpMgHPAm/q9frOGOn5wut94I3ONAC+hG7zHxUk/D3/5E7MgvyIn6ce6TTlZTTdtOCcV89qclbqOmH+6usaVvDzvNaEc0TQ0fiScR7wbJNz8C4gZIxUCz27vQ/8oTMNgK+bt81/VJDwJR96vX6ey2vcGU7AFm/OyYCpxShkU9P2nBaeda2bQ7Q3vq5hTovtnRVI5lPPLpquZaeEujeNjvTnjJJcevZ0H/ikU+MAmt5MhbRYXtHr9dsMBsMYT9tDRb6mi70a57xQx89p0YIW+KljI3BRZ42k/JDxgabtus5aBmwvRCCQQNCNEU5AgaAbIwyAQNCNEQZAIOjGCAMgEHRjhAEQCLoxwgAIBN0YYQAEgm7M/wfNGemg2ZPAqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "test_x = torch.linspace(0, 1, 51)\n",
    "test_y = torch.sin(test_x * (4 * math.pi))\n",
    "with torch.no_grad():\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "\n",
    "lower, upper = observed_pred.confidence_region()\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax.plot(test_x.detach().cpu().numpy(), test_y.detach().cpu().numpy(), 'k*')\n",
    "ax.plot(test_x.detach().cpu().numpy(), observed_pred.mean.detach().cpu().numpy(), 'b')\n",
    "ax.fill_between(test_x.detach().cpu().numpy(), lower.detach().cpu().numpy(), upper.detach().cpu().numpy(), alpha=0.5)\n",
    "ax.set_ylim([-3, 3])\n",
    "ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
